#include <iostream>
#include <ctime>
#include <cstdlib>
using namespace std;
const int M = 4096;
const int N = 4096;

double m1[M][N], m2[M][N], result[M][N];

void (*mul)(double m1[][N], double m2[][N]);

void init(){
    srand(static_cast<unsigned int>(time(0)));
    for(int i=0; i<M; ++i)
        for(int j=0; j<N; ++j){
            m1[i][j] = rand()%200;
            m2[i][j] = rand()%200;
            result[i][j] = 0;
    }
}

void mul0(double m1[][N], double m2[][N]){
    for(int i=0; i<M; ++i){
        for(int j=0; j<N; ++j){
            for(int k=0; k<N; ++k){
                result[i][j] += m1[i][k]*m2[k][j];
            }
        }
    }
}

// void mul1(double m1[][N], double m2[][N]){
//     for(int i=0; i<M; ++i){
//         for(int k=0; k<N; ++k){
//             for(int j=0; j<N; ++j){
//                 result[i][j] += m1[i][k]*m2[k][j];
//             }
//         }
//     }
// }

// 分块(1024个double)
void mul2(double m1[][N], double m2[][N]){
    for(int i=0; i<4096; ++i){
        for(int j=0; j<i; ++j){
            swap(m2[i][j], m2[j][i]);
        }
    }
    double *a1, *a2;
    a1 = m1[0]; a2 = m2[0];
    for(int i=0; i<4096; ++i){
        for(int j=0; j<4; ++j){
            a1 = m1[i]+j*1024;
            for(int k=0; k<4096; ++k){
                a2 = m2[k]+j*1024;
                for(int t=0; t<1024; ++t){
                    result[i][k] += a1[t]*a2[t];
                }
            }
        }
    }
}

void mul_ikj(double m1[][N], double m2[][N])
{ 
    for(int i=0; i<M; ++i){
        for(int k=0; k<N; ++k){
            for(int j=0; j<N; ++j){
                result[i][j] += m1[i][k]*m2[k][j];
            }
        }
    }
}

void mul_ikj_parallel(double m1[][N], double m2[][N])
{ 
    for(int i=0; i<M; ++i){
        // #pragma omp parallel for
        for(int k=0; k<N; ++k){
            for(int j=0; j<N; ++j){
                result[i][j] += m1[i][k]*m2[k][j];
            }
        }
    }
}

int main(){
    mul = mul0;
    init();
    mul(m1, m2);
}